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Abstract. Numerical schemes for the general relativistic hydrodynamic equations are 
discussed. The use of conservative algorithms based upon the characteristic structure 
of those equations, developed during the last decade building on ideas first applied in 
Newtonian hydrodynamics, provides a robust methodology to obtain stable and ac- 
curate solutions even in the presence of discontinuities. The knowledge of the wave 
structure of the above system is essential in the construction of the so-called linearized 
Riemann solvers, a class of numerical schemes specifically designed to solve nonlinear 
hyperbolic systems of conservation laws. In the last part of the review some astro- 
physical applications of such schemes, using the coupled system of the (characteristic) 
Einstein and hydrodynamic equations, are also briefly presented. 



1 Introduction 

The numerical investigation of the Einstein equations is nowadays an impor- 
tant and fruitful line of research in general relativity. There exist a number of 
mathematical formulations of the gravitational field equations which are on the 
basis of all numerical approaches. The level accomplished in the understanding 
of these equations, both mathematically and numerically, is high enough to al- 
low, in principle, for sound numerical approaches. Nevertheless, apart from some 
remarkable results, such as the discovery of critical phenomena by Choptuik j[| 
or, more recently, the achievement of long-term stable three-dimensional null 
cone evolutions of single black hole spacetimes g| , the numerical investigations 
have only partially succeeded, quite understandably due to the complexity of the 
theory, in providing global numerical solutions of generic spacetimes, especially 
in the presence of curvature singularities. Traditionally, the formulation which 
has received the greatest attention by numerical relativists has been the so-called 
3+1 (ADM) formulation ||J]| (see also the recent review by Friedrich and Rcn- 
dall [|| and references therein). This formulation of the Einstein equations as a 
Cauchy (initial value) problem, along with its multiple variants - hyperbolic (see, 
e.g. and references therein) and conformal reformulations |8||| - is still today 
the workhorse of numerical relativity (for an up-to-date summary of the status 
of numerical approaches in 3+1 see, e.g. [[l0| and references therein.) 

On the other hand, despite being known for about forty years now p]JT^Jl^ | , 
the characteristic formulation of the Einstein equations has however been used 
by much less research groups in numerical relativity (see the recent review by 
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Winicour [14Q. And research programs to integrate numerically the conformal 
equations [15] have started only more recently [ p~6lJT7| , |T8| | . Compared to the 3+1 
formulation, the latter two formulations are best suited to study the confor- 
mal structure of the spacetime - the main topic of the workshop -, and the 
propagation of radiation fields within the spacetime. State-of-the-art numeri- 
cal methodology applied to such two formalisms is comprehensively reviewed in 
the corresponding articles by Bartnik and Lehner (characteristic equations) and 
Frauendiener, Husa and Schmidt (conformal equations) in this volume. Apart 
from some test computations involving matter sources presented in Lehner's ar- 
ticle, those papers are mainly concerned with the integration of the vacuum field 
equations. 

The present contribution to this volume aims, on the other hand, at de- 
scribing the current status of the numerical integration of the hydrodynamic 
equations on curved spacetimes, complementing, to some extent, the contribu- 
tions by the above authors. Even though the description will be rather basic and 
general, I will also present some applications and some recent results obtained 
with the coupled integration of the Einstein and hydrodynamic equations within 
the framework of the characteristic formulation of the gravitational field equa- 
tions. It is worth pointing out that, while initially the characteristic evolution of 
matter was limited to idealized systems such as massless scalar fields, nowadays 
it is mature enough to account for fully hydrodynamical evolutions with perfect 
fluids @@|ip2|gg. 

Admittedly, the motivation to develop the capabilities to perform coupled 
evolutions of the matter fields and the geometry needs really not much of an 
emphasis. In astrophysics general relativity plays a major role in the description 
of compact objects in such diverse scenarios as core collapse supernovae, black 
hole formation, accretion, gamma-ray bursts and coalescing compact binaries. 
With the exception of the coalescence and merging of two black holes - the 
number one problem of nowadays' numerical relativity - all realistic astrophysical 
systems and sources of (detectable) gravitational radiation involve matter. 

The only means to study the time-dependent evolution of fluid flow coupled 
to geometry is through numerical simulations. Some scenarios can be properly 
described in the so-called 'test fluid' approximation, where the self-gravity of 
the fluid is neglected. Nowadays there is a large body of numerical investiga- 
tions in the literature dealing with such hydrodynamical integrations in static 
background spacetimes (see, e.g., references in p5|). Most of these are based on 
the pioneering formulation of the hydrodynamic equations by Wilson p6| and 
use numerical schemes based on finite differences with some amount of artificial 
viscosity. The use of conservative formulations of the equations, and their char- 
acteristic information, in the design of numerical schemes started in more recent 
years f27|. 

On the other hand, time-dependent simulations of self-gravitating flows in 
general relativity, evolving the spacetime dynamically with the Einstein equa- 
tions coupled to a hydrodynamic source, are more scarce. Although there is 
much recent interest in this direction, only the spherically symmetric case has 
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been extensively studied thus far (since the pioneering work of May and White 
in 1966 p8[|). In axisymmetry, fewer attempts have been made, most of them 
devoted to the study of the gravitational collapse of rotating stellar cores and the 
subsequent emission of gravitational radiation |29||30| , ^l[ . The three-dimensional 
efforts arc nowadays mainly focused on the study of the dynamics of relativistic 
stars |3^,[33 34 35|,|36| , with the detailed study of the coalescence of close neutron 
star binaries being the key target [ |37[|3l||3l|f4"o| . These investigations are driven 
by the emerging possibility of detecting gravitational waves in a few years time 
with the different experimental efforts currently underway f4lf| . 

The current article deals with the presentation of the main ideas concern- 
ing a particular kind of the "specialized techniques" (in the language of B| ) 
used to solve nonlinear hyperbolic systems of conservation laws with finite dif- 
ferences. The discussion will be specialized to the general relativistic hydro- 
dynamic equations. These equations - as well as their limiting counterparts in 
Minkowski spacctimc and Newtonian gravity - constitute a nonlinear hyperbolic 
system of conservation laws. For such systems there exist ever increasing sound 
mathematical foundations and accurate numerical methodology, imported from 
Computational Fluid Dynamics. The schemes that will be discussed here are 
the so-called high-resolution shock-capturing schemes (HRSC in the following), 
based upon Riemann solvers and written in conservation form. It is worth notic- 
ing that there are a number of excellent textbooks which deal with this subject 
in great detail, in particular jl3 44 , 45| , ^6] | (see also the contribution by Kreiss in 
this volume). Recent reviews on numerical relativistic hydrodynamics are avail- 
able as well [ 47l 43,E5|1 . The interested reader is addressed to these references for 



a more complete information. 

The article is organized as follows: Section 2 presents the relativistic hydro- 
dynamic equations emphasizing work done on conservative formulations. Such 
formulations are well-adapted to the numerical schemes which are discussed in 
Section 3. Applications of these algorithms are shown in Section 4. Finally, Sec- 
tion 5 closes the article with a short summary. 



2 Relativistic hydrodynamic equations 

The general relativistic hydrodynamic equations consist of the local conservation 
laws of the stress-energy tensor T^ v (the Bianchi identities) and of the matter 
current density J M (the continuity equation): 

V M J" = 0. (2) 

As usual stands for the covariant derivative associated with the four- 
dimensional spacetime metric, <7 M „. The density current is given by J M = pu 11 , 
where v) 1 represents the fluid 4-velocity and p is the rest-mass density in a 
locally inertial reference frame. Greek (Latin) indices run from to 3 (1 to 3) 
and geometrized units G = c = 1 are used in the following. 
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By neglecting non-adiabatic effects such as viscosity or heat transfer, the 
stress-energy tensor of a perfect fluid reads: 

T" u =phu^u v +pg» v , (3) 
where p is the pressure and h is the relativistic specific enthalpy, defined by 

fc=l + e + jjj. (4) 

The quantity e is the specific internal energy. 

After choosing an explicit coordinate system x M = (x°,x l ) the previous con- 
servation equations read: 



g 



—gT^ = -V~9K X T^ , (6) 



where the scalar x° represents a foliation of the spacetime with a family of 
hypersurfaces. Furthermore, \J—g is the volume element associated with the 
4-metric, with g = det(gf AtI/ ), and r^ x are the 4-dimensional Christoffel symbols. 

In addition to the equations of motion (Q) and the continuity equation (||) the 
system must be closed with an equation of state (EoS) relating the pressure with 
some independent thermodynamical quantities, such as the rest- mass density and 
internal energy: 

P = p{p,e)- (7) 
Relativistic hydrodynamic flows were first studied numerically with finite- 



difference schemes and explicit artificial viscosity terms 28 2q|. These terms j49| 
were necessary in order to damp the spurious numerical oscillations associated 
with the presence of discontinuities in the flow solution. Such approaches, albeit 
extensively (and successfully) used in different fields of computational relativistic 
astrophysics (e.g., gravitational collapse, accretion, coalescence of compact bina- 
ries, cosmology), were not able, however, to simulate flows with Lorentz factors 7 
larger than 2 pQ ], for which implicit methods were considered more appropriate. 
More recently, however, the use of artificial viscosity terms in non-grid based 
algorithms such as Smoothed Particle Hydrodynamics, has proven not to have 
such severe limitations fill . 

The study of ultrarelativistic hydrodynamics with explicit finite-difference 
methods underwent a revival with the adoption of conservative formulations 
of the hydrodynamic equations and numerical methodology relying upon the 
hyperbolic nature of such system. Theoretical advances on the mathematical 
character of the relativistic hydrodynamic equations were achieved studying 
the special relativistic limit. In Minkowski spacetime, the hyperbolic charac- 
ter of relativistic (magneto) hydrodynamics was exhaustively studied by Anile 
and collaborators (see |52| and references therein). The so-called high-resolution 
Godunov-type schemes, with low numerical dissipation and oscillation-free rep- 
resentation of discontinuous solutions, based upon either exact or approximate 
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Riemann solvers were extended during the 1990s from classical fluid dynamics 
to relativity [ p7|j5^|^ . |55| . [56| . Nowadays, there exists increasing expertise, both 
theoretical and numerical, to investigate extremely fast flows through accurate 
computer simulations (see |47| for a recent review). 

Traditionally, most of the approaches for numerical integrations of the gen- 
eral relativistic hydrodynamic equations have adopted spacelike foliations of the 
spacetime, within the 3+1 formulation [ p6|j5^ , |33[ | Covariant and conservative 
formulations for ideal fluids, have been presented in ]54| and |plf . From the 
theoretical point of view most of the existing formulations of the relativistic 
hydrodynamics equations are written in terms of quantities measured by an Eu- 
lcrian (hxed) observer. By using Eulerian frame variables (relativistic densities 
of mass, momentum and energy) the equations exhibit a conservation form sim- 
ilar to their nonrelativistic counterparts. In most cases, contrary to Newtonian 
hydrodynamics, fulfilling the (desirable) conservation properties is accompanied 
by a nonlinear recovery process to extract physical (primitive) quantities (such 
as rest-mass density, sound speed, etc) from the conserved quantities forming the 
state vector of the system p^ , ^6pl| (we note that a different approach based 
upon a primitive-variable formulation is given in J57| , p5| ). 

As an example, in the formulation developed by Papadopoulos and Font 
the spatial velocity components of the 4- velocity, u % , together with the rest-frame 
density and internal energy, p and e, are taken as the primitive variables. They 
constitute a vector in a five dimensional space w = (p,u l ,e). The initial value 
problem for equations (|B|) and (^3j) is defined in terms of another vector in the 
same fluid state space, namely the conserved variables, U = (£>, S i , E): 

D = J Q = pu° , (8) 

gi = yOi _ phu u i + pg 0i j ( 9 ) 

E = T 00 = phu°u° + pg 00 . (10) 

With those definitions the hydrodynamic equations can be written as a first- 
order flux-conservative hyperbolic system of conservation laws: 

3(V=gU) dW=gF*) _ 
dx° 

The flux vectors F J and the source terms S are given by: 

P i = (JJ ; T J \!P°) = (pu\phu l u? +pg l \phu°u : > +pg 0: >) (12) 



s = (o,-V^dr^\-V^r^ x ). (13) 

The local characteristic structure of these equations has been presented 
m [|lj. For the other conservative formulations mentioned above such infor- 
mation can be found in Refs. |34]|5(||33| (see also [p8|). The relevance of having 
the wave structure to one's disposal in the development of HRSC schemes will 
become apparent in the following section. 
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3 High-resolution numerical schemes 

The hydro-dynamic equations constitute a nonlinear hyperbolic system of con- 
servation laws. Hence, smooth initial data can turn into discontinuous data (i.e., 
crossing of characteristics in the case of shocks) after a finite time during the 
evolution. Standard finite difference algorithms suffer from important deficien- 
cies when dealing with such systems. Typically, first order accurate schemes are 
too dissipative across discontinuities (excessive smearing) while second order (or 
higher) schemes produce spurious oscillations near discontinuities, which do not 
disappear as the grid is refined. 

Finite difference numerical schemes provide solutions of the discretized ver- 
sion of the original system of partial differential equations. Therefore, conver- 
gence properties under grid refinement must be enforced on such schemes to en- 
sure the correctness of the numerical result (i.e., the global error of the numerical 
solution must vanish as the cell width is diminished). For hyperbolic systems of 
conservation laws, schemes written in conservation form are preferred, since - 
as proven by Lax and Wendroff J59| - if convergence exists, it is to one of the 
so-called weak solutions of the original system of equations. Such weak solutions 
are generalized solutions that satisfy the integral form of the conservation system 
d t XJ + d x F = 0: 

poo r+oo />+oo 

/ / (<5 t U + $ x F(\J))dxdt = - / $(x, 0)U(x, 0)dx, (14) 

JO J — oo J — oo 

for any continuously differentiable test function <P(x, t) with compact support. 
They are classical solutions (continuous and differentiable) in regions where they 
are smooth and have a finite number of discontinuities. 

The class of all weak solutions is too wide in the sense that there is no 
uniqueness for the initial value problem. The numerical method should guarantee 
convergence to the physically admissible solution. This is the vanishing-viscosity 
solution of the "viscous version" of the hyperbolic problem: 

dU 3F(U) d 2 V 

dt + dx ~ V dx 2 ' { b) 

when f] — > 0. Mathematically, this solution is characterized by the so-called 
entropy condition (e.g., the entropy of any fluid element should increase when 
running into a discontinuity). The characterization of the entropy-satisfying so- 
lutions for hyperbolic systems of conservation laws was developed by Lax ]60| . 

The Lax- Wendroff theorem ]5!| does not establish whether the method con- 
verges, for which some form of stability is required. Building upon the Lax equiv- 
alence theorem (see, e.g. |6l[]), the notion of total- variation stability has proven 
very successful, although sound results have only been obtained for (nonlinear) 
scalar conservation laws. The total variation of a numerical solution at time 
t = t n , TV(u"), is defined as: 

TV(u«) = £|< +1 -<|. (16) 

i=0 
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A numerical scheme is TV-stable if TV(u") is bounded at any time and 
for any initial data. Present-day research is focused on the development of high- 
resolution numerical schemes in conservation form satisfying the condition of TV- 
stability, such as the so-called Total Variation Diminishing (TVD) schemes p2| |. 
An additional property that any numerical method should satisfy is monotonic- 
ity (see, e.g. [Q). For scalar conservation laws it has been shown that monotone 
methods are TVD and satisfy a discrete entropy condition. Therefore, they con- 
verge in a non-oscillatory manner to the unique entropy (physical) solution. 

In a conservative scheme the time variation of the mean values of the state 
vector U in a given numerical cell - labelled by index i - is given, in the absence 
of source terms, by the flux differences across the cell interfaces. Mathematically, 
such an algorithm reads: 



where At and Ax stand for the time step and cell width, respectively. 

Historically, in 1959 Godunov |i3) developed the first conservative scheme for 
the classical fluid equations in which the numerical fluxes, F i+ i , at every cell 
interface of the computational grid were computed by exactly solving a family 
of local Riemann problems. Such Riemann problems - the simplest initial value 
problem with discontinuous initial data - arise naturally after the discretization 
procedure of the "continuous" solution by means of piecewise constant approx- 
imations. The Riemann problem is invariant under similarity transformations 
(x, t) — > (ax, at), a > 0. The solution is therefore constant along the characteris- 
tics x/t — const, and, hence, self-similar. It consists of constant states separated 
by rarefaction waves, shocks and contact discontinuities. 

Given a general hyperbolic system, if JJ(x,t) = wji(x/t; U_, U + ) is the weak 
solution of a Riemann problem with initial data U = U_ if x < and U = U + 
otherwise, then the numerical fluxes in Godunov's scheme are given by: 



along the characteristic x/t = 0. 

The exact solution of the Riemann problem was extended to relativistic hy- 
drodynamics in pij (see also p5| ) . Since it involves solving a nonlinear algebraic 
system which can be computationally inefficient, in addition to the approxima- 
tion involved in using piecewise constant data, the use of approximate solutions 
of the Riemann problem were proposed. Hence, if w(x/t; U_, U + ) is such an 
approximation, the Godunov-type schemes are defined pt] , |67t as those in which 
U" +1 are computed as: 



u;< 



■n+l 




(17) 



F i+ i =F(w fl (0;U?,U? +1 )), 



(18) 




-Ax/2 
pAx/2 






(19) 
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and the numerical fluxes in the spacetime computational cell [x t _i, x i+ i] x 
[t n ,t n+1 ] are given by: 

i r +1 

where \J(x i+ i,t) is computed by (approximately) solving a Riemann problem 
at every numerical cell interface U(x i+ i,t) = w(0; C/™,U" +1 ). 

The mathematical and algorithmical developments accomplished in the scalar 
case have been extended to nonlinear hyperbolic systems of conservation laws 
using the so-called local characteristic approach. This technique generalizes the 
original procedure due to Roe p8| by applying the scalar algorithms to any 
of the characteristic equations of the system, after a suitable linearization. At 
each interface i + 1/2 of the computational grid, the Jacobian matrix A of the 
system, A = dF/dXJ is assumed to be constant A i+1 / 2 = A(Uj +1 / 2 ), with 
Ui+1/2 being an average between Ui and Uj+i. The original nonlinear system 
is then rewritten as 9 t U + Ad x \J = 0. The eigenvalues of this matrix are the 
characteristic speeds of the Riemann problem. The approximate Riemann solver 
obtains the exact solution of the linearized system, which can be easily computed 
by solving a system of decoupled, linear characteristic (scalar) equations. The 
properties that the matrix A has to fulfill can be found in [Q for the widely 
used Roe's approximate Riemann solver. 

As an illustrative example, for a second-order upwind, monotone scheme such 
as MUSCL |39|, the expression for the numerical flux function reads: 

= \ (f(U?) + FCU^) - ± r^^jj . (21) 

Index p indicates the dimensions of the system. The quantities w = R _1 U are 
the so-called characteristic variables, R being the matrix whose columns are the 
right-eigenvector expressions of the Jacobian matrix associated with the vector 
of fluxes. Furthermore, A and r stand for the eigenvalues and right-eigenvectors 
of such Jacobian matrix. The "tilde" indicates that all quantities have to be 
computed with respect to the linearized Jacobian matrix A. 

The jumps of the characteristic variables at each cell interface are obtained 
by projecting the jumps of the state-vector variables with the left-eigenvectors 
matrix: 

Aw i+i =R7 + yU^ +1 -Uf). (22) 

The left (L) and right (R) states of the conserved quantities U - at any cell 
interface - are computed from the cell-centered values after a suitable monotone 
reconstruction procedure. The way those variables are obtained determines the 
spatial order of the numerical algorithm and controls, in turn, the local jumps 
at every interface. A wide variety of cell reconstruction procedures is available 
in the literature (see, e.g. [ |6^ , [70| , [7ll ) . 
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The last term in the flux-formula, Eq. (|21|) , represents the "numerical viscos- 
ity" of the conservative scheme. The wave structure of the system is thus used to 
provide the smallest amount of numerical dissipation yielding accurate solutions 
of discontinuities without excessive smearing, avoiding, at the same time, the 
growth of spurious numerical oscillations. 

So far we have only considered onc-dimcnsional systems of conservation laws. 
For multidimensional hyperbolic systems containing source terms a standard pro- 
cedure to apply the above schemes is to use dimensional splitting - computing the 
numerical fluxes along every spatial direction independently -, possibly in com- 
bination with a method of lines. Therefore, for a three-dimensional hyperbolic 
system: 

^ + ^ + ^ + M£U) 

at ox ay oz 

where F, G and H are the fluxes in the x, y and z directions, respectively, the 
dimensional splitting algorithm reads: 

T-m+l _ nAtj2 pAt nAt nAt nAtllr-sn fOA\ 

where Cf,C g and Cu denote the operators associated with the corresponding 
one-dimensional PDEs, i.e. the operators computing the numerical fluxes at ev- 
ery cell interface in a given direction. Furthermore, C s is the operator which 
solves a system of ODEs for the source terms: 

<9U 

^=S(U). (25) 

The state vector U at the final time t n+1 is then computed in consecutive sub- 
steps. On the other hand, in the method of lines the time update of all directions 
- and of the source terms - is done simultaneously. The conservative algorithm 
reads: 



dU idtk (t) F i+ y tk -F 



i— i j,k 



dt 





Ax 








Ay 




1 — H, • k _i 

2 hJi K 2 


Az 



(26) 



where the numerical fluxes are given by 



F i+i,i,fc = F(Ui- qt j t k,V i - q+ x ! j < k,...,lJi + g > j t k), 

**i,3+l,k = G(Ui,j-g,fc, Uij-g+i.fc, Uj^+^fc), 

H ij fc+ 1 = H(Vi t j t k-q, U i:ji k-q+l, ...,Ui,j,k+q), (27) 
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q indicating the stencil chosen to compute these fluxes. Further details about 
multidimensional systems and source terms, with particular emphasis in numer- 
ical schemes for stiff source terms, can be found in [7^ ]. 

We end this section by pointing out that during the last few years most 
of the classical approximate Riemann solvers developed in fluid dynamics have 
successfully been extended to relativistic hydrodynamics. The interested reader 
is referred to [[l8| for a comprehensive description of such solvers in relativistic 
hydrodynamics . 

4 Applications 

We now present some applications of the concepts introduced in the previous 
section. In particular we will show some results concerning the numerical evo- 
lution of the equations of hydrodynamics and the gravitational field within the 
context of the characteristic formulation of General Relativity. But let us start 
first with a demonstration in Minkowski spacetime. 

4.1 Shock tube test 

A standard test to calibrate a hydrodynamics code based on the schemes dis- 
cussed in the previous section is the so-called shock tube problem. This is a 
particular version of a Riemann problem in which the initial states at both sides 
of a discontinuity are at rest. Therefore, the state of the fluid at either side of 
the interface only differs in its thermodynamic quantities such as the density 
and the pressure. 

When the interface is removed, the fluid evolves in such a way that four 
constant states develop. In between each state there can exist one of three el- 
ementary waves: a shock wave, a contact discontinuity and a rarefaction wave. 
As mentioned in the preceding section the exact wave pattern of the (ideal) fluid 
state at any given time was first obtained by Godunov p3[ in Newtonian hy- 
drodynamics. Its generalization to relativistic hydrodynamics was accomplished 
by |6J] (see also |65|). This time-dependent problem provides a simple test of 
the shock-capturing properties of any numerical scheme, its level of difficulty 
depending on the initial data. A comprehensive survey of the behavior of a large 
sample of schemes applied to the shock tube problem is presented in (the in- 
terested reader is also referred to for further information on tests commonly 
used to validate numerical schemes for the hydrodynamic equations). 

The main differences between the solution of relativistic shock tubes and 
their Newtonian counterparts are due to the nonlinear addition of velocities and 
to the Lorentz contraction. The first effect yields a curved profile for the rar- 
efaction fan, as opposed to a linear one in the Newtonian case. The Lorentz 
contraction narrows the shock plateau. These effects, especially the latter, be- 
come particularly noticeable in the ultrarelativistic regime (7> 1). 

For our demonstration we consider a fluid whose initial state is specified by 
Pl = 10 3 and pl = 1 on the left side of the interface and by pa = 10 -2 and 
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Fig. 1. The relativistic shock tube problem at time t = 0.35. Normalized profiles of 
density, pressure and velocity vs. distance for the computed (symbols) and exact (solid 
line) solution. All variables were calculated with a third order scheme on an equidistant 
grid of 400 zones. The initial interface was located at x = 0.5. 



pn . = 1 on the right side. This corresponds to problem 2 of p8[ . The adiabatic 
index of the perfect fluid EoS, p = (T — l)p£, is T = 5/3. An initial jump in 
pressure of five orders of magnitude leads to the formation of a thin and dense 
shell bounded by a leading shock front and a trailing contact discontinuity - a 
blast wave. The post-shock velocity is 0.96c (Lorentz factor 7 w 3.5), while the 
shock speed is 0.986c (7 1=3 6). Resolving the thin shock plateau poses a challenge 
for any numerical scheme. 

Fig. [l] shows the results of the shock tube evolution employing a grid of 400 
zones spanning a domain of unit length. The time of the comparison between the 
numerical and the analytic solution is t = 0.35. We have used a HRSC scheme 
based on the HLLE Riemann solver |^,^] and a parabolic reconstruction pro- 
cedure [|r0[ . The solid lines indicate the exact solution. Correspondingly, the 
symbols represent the numerical approximation for the (scaled) pressure (open 
circles), density (cross signs) and velocity (filled circles). As one can clearly see 
the location and propagation speeds of the different features of the solution are 
accurately captured. The shock plateau can be better resolved by simply increas- 
ing the numerical resolution (see ]7^ , p8[| ). Diffusion-free results obtained with a 
one-dimensional exact Riemann solver are presented in p% . 
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4.2 Gravitational collapse of supermassive stars 

Supermassive black holes (SMBH), with masses on the range 10 6 M Q - 10 9 M Q 
(Mq indicating the mass of the Sun) are commonly found in the center of galax- 
ies |7^ , [75| . Supermassive stars (SMS) have been proposed as possible progenitors 
of SMBH. Such stars can develop a dynamical instability f7(|(77j and undergo 
catastrophic gravitational collapse. 

Recently, the gravitational collapse of SMS was proposed by Fuller and 
Shi |r8[ as a possible model for gamma-ray bursts. The neutrino emission from 
the collapse of a SMS could lead to energy deposition by ^-annihilation v ei n iT ) + 
v e ,(fj.,T) ~ * e ~ + e+ ■ Subsequently, 7-radiation would be produced by cyclotron 
radiation and/or the inverse Compton process. 

Trying to shed some light on the viability of that mechanism, Linke et al |p3| 
have studied numerically the gravitational collapse of spherical SMS using a gen- 
eral relativistic hydrodynamics code. The code is based on the hydrodynamics 



formulation developed by 21 1 . The coupled system of Einstein and fluid equa- 
tions is solved adopting a spacetime foliation with outgoing null hypersurfaces. 
In such framework and in spherical symmetry, the Bondi-Sachs metric |I^ , |l2"[ 
reads: 

ds 2 = — du 2 - 2e 2f3 dudr + r 2 (d6 2 + sin 2 6»d0 2 ). (28) 

The Einstein equations simply reduce to two radial hypersurface equations (ODEs) 
for /3(u,r) and V{u, r): 

f3, r = 2nre 4(3 E, (29) 
V r = e 2li - 8vrr 2 e 4/3 5 r - 4irre 4l3 VE, (30) 

where "," indicates partial differentiation. Correspondingly, the hydrodynamic 
equations are given by expressions (|8])-(|l3]), appropriately particularized to spher- 
ical symmetry, and are solved using HRSC schemes Q. 

In addition, the code developed by |2^] includes a tabulated EoS which ac- 
counts for contributions from radiation, electron-positron pairs and baryonic 
gases, as well as energy losses by thermal neutrino emission. 

A typical simulation of a collapsing SMS is depicted in Figure ^[ The initial 
model corresponds to a 5 x 10 5 M© SMS. The figure shows the radial profiles 
of the evolution of the density, temperature, metric components g uu , g ur , and 
radial velocity. Furthermore, the spacetime diagram at the lower right panel 
shows the local proper time against the location of mass shells enclosing fixed 
fractions of the total mass of the star. The arrow indicates the slope of a lightray 
in Minkowski spacetime. One can see that lightrays are severely delayed close 

1 The reader must be aware of the different meaning of the word "characteristic" in 
the context of the Einstein equations and the hydrodynamic equations: while in the 
former case the characteristic formulation of general relativity refers to a particular 
slicing of the spacetime - by means of null cones - in the latter case the notion of 
the local characteristic approach refers to a numerical procedure seeking to exploit, 
algorithmically, the upwind character of the hydrodynamic equations. 
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Fig. 2. Gravitational collapse of a 5 x 10 Mq SMS. Snapshots of the evolution of the 
density, temperature, metric components g uu , g ur and radial velocity u r . The spacetime 
diagram of the lower right panel shows the location of mass shells (AM = 5 x 1O 4 M ) 
versus local proper time r. The lines intersecting the mass shells are hypersurfaces of 
constant coordinate time u and represent trajectories of outgoing lightrays. A black 
hole forms enclosing ~ 25% of the mass of the SMS. The figure is taken from [p3l. 



to the forming black hole. This black hole forms from the innermost 25% of 
the total stellar mass. The collapse lasts 8 x 10 5 s (~ 9.3 days) and the central 
density increases by a factor of 1.08 x 10 7 . The final configuration becomes highly 
relativistic before the simulation is stopped, with g uu = —119 at the surface of 
the star (g uu — —1.0058 initially). Details about the neutrino emission in such 
an evolution can be found in |^3| . 

This and other simulations have been used in p3| to analyze the possibility 
that collapsing SMS could be progenitors of gamma-ray bursts . The com- 



14 Jose A. Font 



prchensive study performed in |23[ reveals that 99% of the energy produced by 
vv— annihilation is deposited in a spherical layer deep inside the star at a ra- 
dius R v < r < 3R V , where R v indicates the location of the neutrino radiating 
volume. Therefore, only a tiny fraction of the energy is deposited near the sur- 
face of the star where excessive baryon loading could be avoided. As a result, 
ultrarelativistic ejection of matter with Lorentz factors 7> 1, a distinctive fea- 
ture of all gamma-ray burst models, cannot be expected in spherical models. 
The simulations performed in |Q show that the spherical collapse of a SMS 
(M > 5 x 1O 5 M ) does not meet the demands for being a successful central 
engine for a gamma-ray burst. 



4.3 Null cone evolution of relativistic stars 

In [^i) we presented the first results of a program we have recently started to 
study the dynamics of relativistic stars by means of null cone simulations in 
axisymmetry. The final aim is to study the gravitational core collapse problem 
and to compute the associated gravitational radiation (7^,0 . 
For these investigations we use the Bondi-Sachs metric: 

ds 2 = - I — e 2/3 - U 2 r 2 e 2 A du 2 - 2e 2f3 du dr - 2Ur 2 e 2 ">du dO 



r 

+r 2 (e 2 ">d6 2 + e- 2 ^sin 2 9 # 2 ), (31) 

with null coordinate u, radial coordinate r, polar coordinate 9 and the azimuthal 
coordinate which is a Killing coordinate. Using this metric the Einstein equa- 
tions split into hypersurface equations on each light cone (for the fields (3, U 
and V), and one evolution equation (for the field 7, not to be confused with the 
Lorentz factor of the fluid), a wave equation (see also Lehner's article in this 
volume). The hypersurface equations, G\ v — 8irT lL , — 0, read: 

P.r = \r (7,r) 2 - \rRrr, (32) 



\^e 2 ^-^U r \ r = 2r 2 



(sin 2 6 l 7) :I . e 



+ 2 7,r7,e) 



-2r 2 R r 



(33) 



V r = - X -/e 2 ^\U. r f + 

4 2Hsm# 



+e 



203-7) 



1 



(sin 0/3, 



sin t 



+ j.ee +3 cot 07,0 - {P.ef 



-2 lt e(l,e-M-^r 2 e 2 "g AB R AB 



(34) 
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where i? M „ is the Ricci tensor and x A — 
equation for the gravitational field reads: 



Correspondingly, the evolution 



4r(r7) i 



2r<y, r V -r 2 2j, e U + sm9 




sin# 2 



(/3 e ) 2 +sin( 



0,e 



sm ( 



MP + P)uo- (35) 



A remarkable property of the above system of equations is that they form a 
hierarchy: knowing 7 on the first null hypersurface allows one to radially integrate 
the corresponding equations to determine 0, U, V and 7 iU (in that order) on 
that hypersurface H]. 

In the code developed by p4| the numerical implementation of the Einstein 
equations closely follows that of |30|| . The same marching algorithms are em- 
ployed with additional source terms arising from the presence of matter fields. 
Since the code uses spherical coordinates, special care is taken with the nu- 
merical treatment of the coordinate singularities at the origin and at the polar 
axis. In order to impose boundary conditions at the origin the assumption that 
t = u + r, x — r sin 9 cos <j), y = r sin 6 sin (f> and z = r cos 9 form a local Fermi sys- 
tem at r = is enforced. This implies a fall-off behavior given by V = r + 0(r 3 ), 



(3 = 0(H), U = 0(r) and 7 = 0(H) 



Regularity on the axis requires that 



XJ I sin# y 7/ sin 9 are continuous functions at 6 = 0,7r. The code also uses a 
new polar coordinate y = — cos 9. In order to keep the freedom of working with 
numerical grids which only cover the star without its vacuum exterior, the radial 
coordinate used in |SC| was generalized: Starting from an equidistant radial coor- 
dinate x £ [0, 1], the code of |2^] allows for a general coordinate transformation 
of the form r = r(x), so that either compactified - with future null infinity being 
the outermost radial grid point - or non-compactified grids can be used. This 
will allow for the unambiguous computation of the gravitational radiation at T + 
in our planned gravitational core collapse simulations, by simply reading off the 
news function at the outermost radial grid point. The hydrodynamic equations 
are formulated as in (2^] and solved using HRSC schemes. 

As a simplified model for a self-gravitating relativistic star |24| have con- 
sidered the spherically symmetric solution of the general relativistic hydrostatic 
equation, the so-called Tolman-Oppcnhcimer-Volkoff equation, with a polytropic 
EoS p = Kp r . Equilibrium models were used to check the long-term stability of 
the code and its convergence properties. The code has proven to be stable for 
evolution times much longer than the characteristic light-crossing times of the 
different models considered. Spherically symmetric simulations have also shown 
the expected second order accuracy of the code. 

Following |36| we have also checked the code on a dynamical evolution of 
an unstable spherical initial model. In such a model the sign of the truncation 
error of the numerical scheme controls the fate of the evolution. In the code 
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Fig. 3. Evolution of the central rest-mass density during the migration of an unstable 
relativistic star (JV = 1,K = 100, M = 1.447Af , p c = 8.0 x 10~ 3 ;G = c = M Q = 1) 
to a stable model with the same rest-mass. The central density of the (final) stable 
configuration is p c — 1.35 x 10~ 3 . The evolution shows the expected behavior. Since we 
are using a polytropic EoS, the amplitude of the oscillations is essentially undamped 
for the evolution times shown. 



this sign is such that the unstable star "migrates" to the stable branch of the 
sequence of equilibrium models. In such a situation, the rest-mass of the star has 
to be conserved throughout the migration. Despite the fact that this mechanism 
cannot occur in nature - unstable stars can only collapse to more compact con- 
figurations - and as such it is an academic problem, it represents, nevertheless, 
an important test of the accuracy and self-consistency of the code in a highly 
dynamical situation. 

As in ||1 we have constructed a N = 1 (T = 1 + 1/N = 2), K = 100 
polytropic star with mass M = 1.447 M and central rest-mass density p c = 
8.0 x 10~ 3 (in units in which G = c = M Q = 1). Fig. || shows the evolution of 
the central density up to a final time of u = 1500. On a very short dynamical 
timescale the star rapidly expands and its central rest-mass density drops well 
below its initial value, less than p c — 1.35 x 10~ 3 , the central rest-mass density of 
the stable model of the same rest-mass. During the rapid decrease of the central 
density, the star acquires a large radial momentum. The star then enters a phase 
of large amplitude radial oscillations around the stable equilibrium model. As 
Fig. U shows the code is able to accurately recover (asymptotically) the expected 
values of the stable model. Furthermore, its evolution is completely similar to 
that obtained with an independent fully three-dimensional code in Cartesian 
coordinates |36| . 

The evolution shown in Fig. |^ allows to study large amplitude oscillations 
of relativistic stars, which cannot be treated accurately by linear perturbation 



Local characteristic algorithms for relativistic hydrodynamics 



17 



1e-05 



8e-06 



6e-06 



4e-06 



2e-06 




0.0025 0.003 0.0035 0.004 0.0045 0.005 
radial resolution (1/grid zones) 



Fig. 4. Convergence properties of a global energy conservation test. The difference 
between the initial and final Bondi mass converges in an almost second order way to 
the value given by the energy which is carried away by gravitational waves. The figure 
is taken from |24|. 



theory. These oscillations could occur after a supernova core-collapse or after 
an accretion-induced collapse of a white dwarf. 

To end this section we briefly discuss a global energy conservation test of 
the axisymmetric characteristic code which was presented in p4j| . For that pur- 
pose Siebel et al |24| used a strong ingoing gravitational wave to perturb an 
equilibrium relativistic star: 



1- 



2- = 0.05 e- 2 ^ 4 ) 2 e-^" 



sin 6 2 



(36) 



The (nonlinear) initial pulse induces large velocities in the fluid of the star 
which give rise to "strong" outgoing gravitational waves, i.e. with an energy 
larger than the numerical errors involved in the calculation of the Bondi mass 
for a given resolution and integration time. If M is the Bondi mass and P is the 
total energy radiated away by gravitational waves at future null infinity X + , the 
convergence of the quantity 



ec := Ml 



M 



u—u^ >0 



p\ 



[0,«, 



(37) 



to zero represents a very severe global test of the numerical code. Satisfactory 
convergence results under grid refinement are shown in Fig. ^. 



5 Summary 



The article has dealt with presenting some concepts and applications in rela- 
tivistic astrophysics of a particular class of finite difference numerical schemes 
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based on Riemann solvers, specifically designed for nonlinear hyperbolic systems 
of conservation laws. 

Such schemes have been discussed in the context of the general relativis- 
tic hydrodynamic equations. Nevertheless, the algorithms presented are gen- 
eral enough to be applicable to other hyperbolic systems such as the Einstein 
equations (when appropriately formulated, as e.g. in Friedrich's conformal ap- 
proach JL5| ; see also While this may not be strictly necessary for vacuum 
spacetimes, it may become relevant when dealing with nonvanishing stress- 
energy tensors. 

The use of conservative algorithms based upon the characteristic structure 
of the hydrodynamic equations, developed during the last decade building on 
ideas first applied in Newtonian hydrodynamics, provides a robust methodology 
to obtain stable and accurate solutions even in the presence of discontinuities. 
This has become apparent since the early 1990s 

The knowledge of the wave structure of the equations is the essential build- 
ing block in the construction of the so-called linearized Riemann solvers. The 
increasing use of these solvers in relativistic hydrodynamics has proved suc- 
cessful in handling complex flows, with high Lorentz factors and strong shocks, 
superseding more traditional methods based on artificial viscosity |47j . 

In the last part of the article we have discussed some astrophysical applica- 
tions of such schemes, using the coupled system of the (characteristic) Einstein 
and hydrodynamic equations. Examples involving the gravitational collapse of 
supermassive stars and the evolution of relativistic compact stars have been 
presented. 
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